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Abstract 

We consider a stochastic SIS infection model for a population partitioned into m 
households assuming random mixing. We solve the model in the limit m — > oo by 
using the self-consistent field method of statistical physics. We derive a number of 
explicit results, and give numerical illustrations. We then do numerical simulations 
of the model for finite m and without random mixing. We find in many of these 
cases that the self-consistent field method is a very good approximation. 
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1 Introduction 



In this paper we will give a treatment of an SIS infection model in a population 
partioned into households. Our main emphasis will be the use of a method 
borrowed from statistical physics, the self-consistent field (SCF) technique. 

The susceptible-infected-susceptible (SIS) model is one of the basic themes 
in mathematical epidemiology [1,2]. In this idealization of the spread and 
persistance of an infection, N individuals are initially either susceptible, S, or 
infected, /. Each one of the infected is coupled to a certain number of other 
agents. Each / infects all of the S with which it has contact with a rate B, 
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and each I recovers with rate 7, and immediately becomes susceptible again. 
Schematically, we can write: 



S + I 



21 



I 



S 



(1) 



Clearly, at long times the system might achieve an equilibrium where the 
number of infected fluctuates around a mean. We characterize this state by 
denoting the mean fraction of population which is infected by /, i.e., 



where <> t denotes a time average. We refer to this quantity as the endemic 
level. 

The model is defined by the two parameters, B, 7, and by the contact structure 
of the population. Two cases of the latter have been studied in some detail. 
In the first, each individual is equally coupled to all the others in a large 
population. This random mixing model is rather easy to treat, and a large 
number of results are known [2]. For example, in the deterministic limit, I 
obeys: 



where we have redefined B = b/N. It is easy to show that for b > b c = 7, and 
for any intial condition with I > 0, the epidemic goes to the endemic level, 
/ = 1 — 7/6. For b < b c the infection dies out. 

In the second well-studied case the individuals are arranged in some geometric 
way and couple only to their nearest neighbors. This is called the contact 
process [3]. In the last section of this paper we give numerical results on a 
model which interpolates between the two cases. 

Ball [4] introduced the SIS household model, which puts more structure into 
the random mixing case by having two levels of mixing. In this idealization 
the population, N, is partitioned into m households with n members. The 
infection rate within the household, w, is different from that with individuals 
outside, B. This followed work on SIR models with similar structure [5]. Ball's 
results are for the case m — > 00 with fixed n, that is, a large population of 
small households, which is clearly of practical interest. 

In this paper we extend Ball's work from a different point of view inspired by 
techniques in statistical physics. These will be explained in the next section. 
Briefly, the idea is that for random mixing and a large population, we can 



f=<I>t /N, 



(2) 



dI/dt = bIS/N -7/ 
S = N-I 



(3) 
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consider that each household is acted on by an average force of infection from 
outside. This is quite similar to methods used in the theory of magnetism 
[6], inspired by the work of P. Weiss. Weiss introduced the 'molecular field', 
which is proportional to the mean orientation of the neighbors of a given 
magnetic ion. The orientation of the ion itself is then determined, and the 
unknown molecular field is gotten by demanding that the neighbors have the 
same orientation. (Note that this step assumes that there is no correlation 
between neighbors.) This step is known as self-consistency, hence the name, 
self-consistent field. As we will see, / plays the role of the self-consistent field. 
The technique itself is called self-consistent field (or sometimes mean-field) 
theory. 

Some of our results were already obtained in [4], e.g., the expression for b c for 
households of size n, Eq. (27), below. Other explicit results are new. We think 
that one of the main values of this paper is that this point of view allows us 
to look at the solutions in a new way, and gives considerable insight. 

In the last section we continue the analogy to magnetic systems by doing 
numerical simulations for finite to, and for systems where the random-mixing 
assumption is not valid. We find, as in magnetism, that the SCF assumption 
is remarkably accurate in situations where we would expect it to fail. 



2 Stochastic formulation and steady-state solutions 

Consider an SIS epidemic in a population partitioned into to households each 
consisting of n individuals. We suppose that within a household the infection 
rate per contact is w and between different households B. We write B = 
b/[(m — l)n]. (In the notation of Ball [4], \ w /n = w, A& = b ). We suppose that 
b remains finite in the limit to — > oo. The recovery rate is 7, and a recovered 
individual can be reinfected at once. 

2.1 Self- consistent field equations 

We begin intuitively by considering the case n = 2 and try to formulate the 
Chapman-Kolmogorov forward equations for the probability, Pk, to have k 
infected in a typical household: 

7P1 - 2bfP 

2bfPi-( 1 + w + bf)P 1 + 2 1 Pi 



dPo 
dt 
dPi 

~dt 



3 



dP 2 



f = (w + bf)P(- 
f = (P 1 + 2P 2 )/2 



-2 7 P 2 



(4) 
(5) 



The terms in b are the mean force of infection between households, and / has 
the interpretation of the fraction of a typical household that is infected. The 
point is that the only communication between households is by infection, and 
for a large number of such contacts, only the average number of infected in the 
other households is important. In effect, we have assumed that the households 
are independent. Note that J2k Pk = 1, and < / < I. 

The strategy for solving these equations (in the steady state) is to find the Pk 
with / as a parameter, and solve for / self-consistently, using Eq. (5), at the 
end of the calculation. This is a heuristically appealing method. For magnetic 
models, the approximation is exact in the case of infinite-range interactions (in 
our case, all households interact equally with all others). Even for short-range 
interactions (i.e., households interacting with a 'neighborhood') the approxi- 
mation is often surprisingly useful. 

In fact, a rigorous proof of the validity of Eqs. (4, 5) is available for this case 
as well. This is given by Ball [4] who uses the work of Kurtz [7] to prove the 
following: let X^ be the number of households in the population with exactly 
k infected. Define Pk = lim,^,^ xj^/m. Then P k satisfies Eqs. (4, 5) for any 
initial condition with finite density, i.e., with P k (0),k > not all equal to 0. In 
the last section we will investigate numerically how good this approximation 
is for finite m and for finite range of interaction. 

We note, for later use, that we can give another interpretation to Eq. (4). 
Suppose we consider a single household with 2 individuals subject to a fixed 
external force of infection. Then we can think of P k as the probability to have 
k infected. From standard results on stochastic processes, there is a unique 
steady-state solution to these equations for any positive /, b, w. Further, if 
/ = the only steady-state solution is P = 1, Pi — Pi — 0; that is, the 
infection becomes extinct with probability unity. 

In the general case the time evolution of the vector of probabilities P = 
(Po, Pi, P n ) T ( T denotes a column vector) to have exactly k infected in the 
household of n, is given by: 



dP 

~dt 



M(/)P 



(6) 



where M(/) is a square tridiagonal matrix with nonzero elements 



m k ,k-i — (n — k + i)[(k — i)w + fb] 



(7) 
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Tfik.k = — [kj + (n — k)kw + (n — k)fb] 
m k;k+1 = (A; + 1) 7 

(0 < k < n). The equations again contain /, the mean level of infection, as a 
parameter. 

1 n 

f = -Y / kP k (8) 

n k=0 



For the reasons given above, these equations are valid in the limit m — > oo. 

We will be interested in the stationary solutions of Eqs. (6, 8). It is clear by 
inspection that there is always a trivial steady-state solution / = 0, Pq = 
1, P k — 0, k > corresponding to extinction of the infection. We will find that 
above a certain threshold there is another solution, and we will show that 
when this occurs the trivial solution is unstable. 



2.2 Steady-state solutions 



In the case n = 2 it is a matter of simple algebra to find the steady-state 
solution: 



Pl 2fb -< 



7 2 + 2/6 7 + (fb) 2 + fbw 

(fb) 2 + fbw 
^ + 2fb 1 + (fb) 2 + fbw 



(9) 



Using Eq. (5) we get a cubic polynomial equation to solve for the endemic 
level, /. After removing the trivial solution we have: 

(fb) 2 + /(2&7 + wb - b 2 ) + ( 7 2 - 67 - wb) = (10) 



The larger root of this equation is: 

(6 - w - 27) + J \b + w) 2 + 4w7 
/ = - (11) 



This expression is acceptable, i.e., positive, provided: 

b( ri + w)h 2 >\ (12) 



5 



Thus Eq. (12) serves as a threshold condition for the existence of a non-trivial 
endemic level. 

For n = 3 it is not difficult to find a cubic equation for / analogous to Eq. 
(10) (in this expression we set 7 = 1, which sets the unit of time): 

fb 3 - f 2 (b 3 - 3wb 2 - 3b 2 ) - f(b(2b + 3bw) - 2bw 2 -3b- 3bw) 

-b(l + 2w + 2w 2 ) + l=g(f) = (13) 

By inspection g(l) > 0. If g(0) = —b(l + 2w + 2w 2 ) + 1 < then there is at 
least one root in the interval [0, 1]. (In fact, numerics shows that there are also 
two negative roots). Thus if: 

b(l + 2w + 2w 2 ) > 1 (14) 



we have a non-trivial solution, and otherwise all the roots of Eq. (13) are all 
negative, and only the trivial solution is acceptable. 

In the general case we look for the stationary solution of Eq. (6) by setting 
dP/dt = 0. We get a homogeneous system of equations, from which the en- 
demic level has to be established self-consistently. Recall that J2k Pk — 1, so 
that we can reduce the number of equations by one by eliminating P . Denote 
the vector of length n, p = (Pi, ...,P n ) T . The resulting system of equations 
for the stationary state is inhomogeneous. It reads 



S(/)p = -n/6u 

u=(l,0,0,...) T 



(15) 



The matrix S is given by: 



m i,i — nfb mi t 2 — nfb —nfb —nfb... —nfb 
m 2 i m 22 m 2 3 ... 

... m fejfc _i m fejfc m Kk+ i 0... 



V 







m„ )n _i m nn j 



(16) 



The matrix S(/) is linear in /: 
S(/)=W + /6B 



(17) 
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where both matrices, W and B are /-independent. The structure of W is: 



-j-(n-l)w 2 7 ... ^ 

(n-l)w -27-2(n-2)iu 37 

2(n-2)w -37-3(n-3)w 47 

3(n-3)w 

... (n — l)w —717 J 



(18) 



W depends only on the parameters 7 and w, characterizing the infection 
within a household. Note that it is nondegenerate: detW = (— 7) n n!. To 
see this, add to each row the sum of all the rows below. The result is a 
lower triangular matrix with diagonal elements in the n-th row equal to —717. 
Moreover, B is a matrix with integer components: 



B = 



— (2n — 1) — n —n —n 

(n-1) -(n-2) 
(n-2) -(ra-3) 



\ 



1 



(19) 



We give the explicit forms of the matrices W and B for n 
setting 7 = 1: For n = 2 one has 



W 



— 1 — w 



w 



-2 



and B 




2 and n — 3, 



(20) 



For n = 3 we get 
/ 

W = 



-l-2w2 
2w -2 - 2w 3 

2w -3 



^ -5 -3 -3^ 



and B 



2 -10 
1 



(21) 



From the remarks above, it is clear that there is a non-zero solution, p, to Eq. 
(15) for any / > 0, so that S -1 exists. Thus: 



p = -n/6S _1 u 



f22) 
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From this equation a closed form for the endemic level / follows. Multiplying 
both sides of Eq.(22) by the vector rix = (1,2, ...,n) and noting that / = 
n _1 n!p we get: 

/ = -fbmS- l u (23) 



which is a closed algebraic equation for / which always has the trivial solution 
/ = 0. Factoring out this out, we get an equation for the nontrivial solution: 

- 6niS _1 u = 1. (24) 



As in Eqs. (12, 14) there is a condition on the parameters for the solution to 
be positive. For given w we define a value b c such that for b > b c we have such 
an /. This can be found from Eq.(24) by noting that at b c the second solution 
crosses / = 0. At that point, S = W, and we get: 



b c = -1/ mW 



u 



(25) 



To get an explicit expression for b c consider the vector v = W 1 u, whose 
elements are v fc = [W -1 ^ k . As we will show, 



(n-iy. fe _i 



(26) 



so that the threshold is given by 



S.-1 



^ (n-l)\ fc _! 

-mv = > 7 -4tw • 



(27) 



To prove Eq.(26) we note that the solution to Wv = u is unique so that we 
need only insert Eq. (26) into Wv = u . The first element is 

- (-1 - (n - l)w) - (n - l)w = 1, (28) 



while the k-th element (k > 1) is: 



(k - l)(n + 1 - k)w 



(n-l)\ 



(n + 1 - k)\(k - 1) 



w 



k-2 



in - IV 

+ (_A; _ k (n - k)w)± 7TTT-U' + 

(n — ky.k 
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+(k + l) 

(n-l)\ 

+ 
0. 



(n-l)\ 



n-k- l)\(k + 1 
1 1 



-w 



(29) 



(n — k)\ (n — k)\ 
1 



(n-k-l)\ (n-k-l)\ 



w^ 1 + 



w 



Once we obtain p we can find other moments of the distribution by multiplying 
Eq.(22) by n q = (1, 2 q , 3 q , n q ). We present numerical results on n 2 p below. 

We now discuss the behavior of the solution in the vicinity of the transition 
point. Eq.(22) can be rewritten: 



p = - n fb (i + fbW^B) 1 W- 1 u ! 



(30) 



Then we formally expand in a geometric series: 
p = - nfb (l - fbW^B + (fbf (W _1 B) : 



... ) W- X u. 



(31) 



Multiplying both sides of Eq.(31) by n 1 ni we get the following equation for 
/: 



/ = -fb (a - fba, + {fbfa 2 + ... + (-l) m (fb) m a m + ...) 



(32) 



where 

a m = ni (W-^)" 1 ^^-^ 



(33) 



from which / can be obtained with any necessary accuracy. Note that accord- 
ing to Eq.(25) the coefficient 



a = niW x u = — l/b c . 



(34) 



Eq.(32) is especially useful in the vicinity of transition, where / is small. 
Both the point of transition and the endemic level near the threshold can be 
determined. 

Eq.(32) always possesses a trivial solution / = 0. An equation for a nontrivial 
one then reads: 



G(f) = fb' ai - f b A a 2 + ... + {-l) m+1 f m b m+1 a m + ... = 1 + ba . (35) 
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The formal solution of this equation involves the inverse of G: 
f = G-\l + ba ) = G- 1 (l-b/b c ). 



(36) 



where Eq (34) has been used. Using the relationship between the Taylor series 
of a function and its inverse we get, in the vicinity of b c : 

f = -(3 lT + (3 2 T 2 -... (37) 



where r = {b — b c )/b c and the coefficients are: 



^ 6 2 ai(w)' 

$2= — i^L "T j (38) 



^ 3 = m S ( 2 °2( w ) - a 1 (w)a 3 (w)) , 
fl-^ [ tv J 



In the vicinity of the transition one can take b = b c (w); the behavior of / is 
dominated by the linear term: 

/ = -P ± r = -j^T^r (39) 
b 2 ai(w) 



Explicit calculations show that a\ < 0. 



2.3 Stability of the disease-free equilibrium 



Restoring the time- dependence in Eq.(15) we obtain an equation for the time 
evolution of the reduced probability vector p: 

Ap = (w + /6B)p + nfbM. (40) 

(XL 



Since p depends on /, Eq.(40) is a nonlinear equation describing the evolu- 
tion of p. In the previous section we have shown that this equation has two 
equilibria provided b > b c , where b c is defined in Eq. (27). We will now show 
that under the same conditions the trivial solution, / = 0, p = is unstable. 
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In the vicinity of the trivial equilibrium, we can omit the term /6B in brackets 
in Eq. (40). The linearized equation reads 



d 
dl 



p = Wp + 6(nip)u = Tp 



(41) 



where 

T = W + 6N, 



(42) 



The matrix N has rix as its first row and all other elements zero. 

The determinant of T is linear in b. To see this, note that the explicit form of 
T is: 



-7-(n-l)iu + & 27 + 26 36 
(n-l)w -27-2(n-2)iu 37 

2(n-2)w -37 - 3(n - 3)w ... 



\ 



in — l)w 



(43) 



-717 J 



If we expand the determinant about the first row, each term contains b only 
once. It can thus be written det T = det W + cb, with c a constant. However, 
the determinant of T is the product of its eigenvalues. Since for b = (separate 
households) the trivial solution is stable, all of the eigenvalues of W are neg- 
ative; det W is nonzero, and its sign is positive for n even and negative for n 
odd. From the linearity in b we see that det T changes sign at b = —c" 1 det W 
which corresponds to the appearance of a positive eigenvalue. 



To find the b for which det T = write Eq.(42) as: 



TW 



I + &NW~ 



(44) 



The determinant of the matrix in the l.h.s. of this expression is 



det T/ det W = [(- 7 )" n /n! 



detT 



(45) 
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Fig. 1. The endemic level as a function of b for n = 3 and various w, 7 = 1 

and vanishes when det T = 0. But, since N has non-zero elements only in the 
first row, the matrix on the r.h.s. of Eq. (44) is: 



1 + 691 bq 2 bq 3 

10... 

10 



bq, 



\ 



(46) 



where is the i-th element of the vector niW -1 . By inspection, the determi- 
nant of this matrix is 1 + gi = l + 6n 1 W~ 1 u. It vanishes at 6 = —1/ (niW _1 u), 
i.e. exactly at 6 = 6 C , see Eq.(25). Therefore, for 6 > 6 C a positive eigenvalue 
appears, and the trivial solution loses its linear stability. 



2.4 Numerical results 



It is quite simple to work out the dependence of / on the parameters w, 6, n. 
We wrote a Matlab program which solved Eq. (24). In Figure (1) we show / 
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0.2 0.4 0.6 0.8 



Fig. 2. The standard deviation of the number of infected in a household, as a function 
of b for n = 3 for several w's, 7 = 1 



- 0.5 




Fig. 3. The endemic level as a function of b for w = 0.2 and n = 1,2, ..8, 7=1 
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0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 



Fig. 4. The standard deviation of the number of infected in a household, as a function 
of b for w = 0.2 and n = 1, 2, 8, 7=1 

as a function of 6 for several u>'s for n = 3. In this and the following figures 
we have taken 7 = 1. For other values of 7 one can use the obvious scaling 
relation (w,b, 7) — > (10/7,6/7,1). 

Recall that / = rixp can be interpretated as < i > /n where i is the number 
infected in a household, and <> is the average over the distribution P. If 
we calculate n 2 p we get a measure of the fluctuation of i in a household. In 
Figure (2) we show the standard deviation, a = [< i 2 > — < i > 2 ] 1 / 2 = 
[n 2 p — {nf) 2 ] l l 2 as a function of b for several u>'s, for n = 3. In a similar way, 
in Figures (3, 4) we show / and a for w = 0.2 and n = 1, 2, 8. 



3 Simulations: finite systems and finite range 

The considerations in the previous sections apply in the limit m — > 00. For 
any finite m the only steady-state solution possible is the trivial one, / = 0; 
with probability unity the infection will become extinct. However, the time to 
extinction may be very long, 0(e m ) if the system is above b c . See, for example, 
[8,9,10] for the case n = 1. For times short compared to the extinction time the 
system will be in a quasi-stationary state which resembles the infinite system. 
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Fig. 5. The coupling scheme used in the simulations. This is the case n = 3, m = 8, 
L = 2. The white bonds represent contacts within a household. All of the in- 
ter-household contacts for one agent are shown in black 

We have done preliminary investigations of this and related behavior using 
computer simulations. We represent the population as follows: imagine a pop- 
ulation of m households each containing n agents arranged on a ring. All of the 
members of the same household can infect one another with probability/unit 
time w. The communication with other households is variable range: we allow 
each household to infect all of the other households in a range L on either 
side. See Figure (5). That is, the number of outside contacts is q = 2Ln, and 
the infection probability is b/q. The equivalent of the random mixing case of 
the previous sections is 2L — m — 1; in this case all agents are coupled to each 
other. The case n — 1, L — 1 is the well-known contact process [3]. 

Our computer algorithm is quite simple: we suppose that 7, w,b are all less 
than unity. We start by marking a certain fraction of agents as infected, which 
establishes 1(0) where I is the total number infected in the population. For 
each computer time step we choose one agent, in household k, at random 
among the infected, and find all the other agents that are coupled to the one 
in question. These will be in households with index k — L, k — L + l, ...,k + L, 
where the indicies are to be interpreted with periodic boundary conditions. 
With probability w or b/q (depending whether the other agent is in k or 
outside) the other agent is infected. Then another agent is chosen at random 
among the infected, and it is allowed to recover with probability 7. Each 
computer time step corresponds to an increment of 'clock time', dt = 1/1. 
A typical result for I(t) is shown in Figure(6). It is worth noting that this 
algorithm is quite fast - it is not difficult to treat as many as 3000 agents on 
a simple workstation. 

In this section we will give a few results of the simulations. We have found 
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Fig. 6. Simulation results for the number of infected as a function of time compared 
to the SCF results, n = 3, m = 501, L = 150, w/j = 0.2,5/7 = 1. The agreement 
for the average is essentially perfect, but for the fluctuations (see text) the SCF 
method is an underestimate. 

that the SCF is an excellent representation of the behavior in large systems 
with long-range interactions, but that for smaller systems and shorter ranges 
the approximation becomes poorer, but still surprisingly good. A systematic 
investigation of these effects will be reserved for a future publication: here we 
give a few preliminary results. 



3.1 Finite size effects 



For large enough m the simulation results agree very well with the SCF results. 
In Figure(6) the average of the time series agrees almost perfectly with the 
results from the previous sections: the numerical average of the simulation 
data differs from the SCF by less than the width of the line. 

In addition we can try to understand the fluctuations of the time series by the 
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following reasoning: if we imagine that the different households fluctuate more 
or less independently (this is the assumption under which the SCF method is 
valid), then the standard deviation of the total should be a tot = a^/m. As we 
see in the figure, this is a good approximation. However, we have observed that 
there is a large finite-size effect for this statistic. For to = 301 the simulation 
result for a tot is almost twice the prediction. 

Another way to look at the standard deviation is to simply record the time 
series of the number infected in a given household, and directly find a. We 
did this by averaging over 10 randomly chosen households in the system. 
We find very good agreement with the analytic treatment: for example, for 
w = 0.2, 6 — l.,n — 3, to = 501 we find a = 0.82. The analytic result is 0.84, 
well within our estimated error. 

As we reduce to we find a surprising result: within the numerical accuracy of 
our simulation we can find no change in /. Rather, at a certain small size, 
to = 50 for the parameters in the last paragraph, the infection dies before we 
can collect adequate statistics. 

3.2 Finite range effects 

For the case of small L the SCF ceases to be valid. In magnetic systems, 
locality of the interaction is well known to increase fluctuations and to change 
the threshold. We find related effects here, but the range needs to be very 
small indeed for the SCF result to be very inaccurate. 

In Figure (7) we show simulation results for f(L), the endemic level as a 
function of the range of the interaction between households, L. We show three 
different sets of parameters for n = 3 and compared with the SCF (L — > oo). 
The dependence is remarkable: for the parameter values (w — 1., 6 — 1.), finite 
range has hardly any effect. The other two sets do show a decrease of /, and an 
apparent shift in the threshold. For w = 0.2 b c = 0.67 in the SCF, well below 
6=1. However, for L < 5, the system is below threshold, and the epidemic 
dies quickly. A similar effect occurs for w — 1., 6 = 0.3, which is below 6 C for 
L < 10. (For large L, 6 C = 0.2 for these parameters.) Another way to put this 
is that there is a characteristic length L c ps 10 for these two sets of parameters. 

These results are reminiscent of effects that occur in phase transitions in sta- 
tistical physics. In fact, the contact process is fruitfully viewed as an example 
of a continuous non-equilibrium phase transition [11]. Near the threshold (i.e. 
the critical point) of a continuous phase transition, the system develops a di- 
verging correlation length. For total system sizes greater than this length finite 
size effects are small. Thus we might guess that the correlation lengths for two 
of the sets of parameters could be of order L c . We have not investigated this 
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Fig. 7. The endemic level, /, as a function of the range of the infection contacts 
between households. The horizontal lines are the SCF results for large L. In all cases 
n = 3, m = 501 so that the maximum possible L = 250 

matter in detail, but it would repay consideration. 



4 Summary 



In this paper we have given a number of explicit results for the SIS household 
model. Our point of view was to emphasize the SCF method as a guide to 
solution and insight. The numerical results give rise to a number of questions 
that are probably not accessible to rigorous proof, but for which numerical 
methods could be pursued. For example, we have speculated that there are 
critical fluctuations [6,11] in this system. That would mean that near threshold 
there would be large correlations in adjacent households for a finite range 
system, a result that could be of considerable interest. 

There is also a practical aspect to this method. SCF theory is quite flexible, 
and allows considerable complication to be added to the model without sub- 
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stantially increasing the complexity of the solution. Suppose our household 
contains two types of members, for example, children who go to school and 
contract infections, and adults who do not. This would require a few changes: 
for a 3 person household with two adults and a child, the analog of the vector 
P would be of length 6 instead of 4, and the complete solution would be no 
worse than the numerical inversion of a 5 x 5 matrix. 
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